1 Introduction

In this script, we go through all the pre-registered proposed analyses. As a reminder, the main research questions where as follows:

  1. To what extent does infants#’ preference for IDS as measured in a laboratory setting predict their vocabulary at 18 and 24 months?
  2. Does the relation between IDS preference and vocabulary size change over development?
  3. Are there systematic differences in the strength of this relationship across the language communities in our sample?

We first present the main “sample theory based” analyses (also known as frequentist), separately on the North American and UK samples in parallel to answer our first two research questions, then on the total dataset to answer our third research question. We then provide additional Bayesian statistics where a null effect was found, as specified in the pre-registration.

# Library imports, general settings ==============
library(tidyverse); library(egg)
library(knitr)
library(lme4); library(lmerTest); library(simr)
# As in our discussion with Mike, we will use lmerTest for calculating p values
# in the mixed-effects models (noted as a deviation) 
library(brms); library(rstan)
library(future)
plan(multicore, workers = 8) # Swap to multiprocess if running from Rstudio
library(lattice)
library(effects)
library(sjPlot)
library(robustlmm)
library(car)

# Load model comparison functions
source("helper/lrtests.R")

# Deal with package priority issues
select <- dplyr::select

theme_set(theme_bw(base_size = 10))
options("future" = T)
#knitr::opts_chunk$set(cache = TRUE)

print(sessionInfo()) #listing all info about R and packages info
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] furrr_0.2.1          future.apply_1.7.0   bridgesampling_1.0-0
##  [4] car_3.0-10           robustlmm_2.4-4      sjPlot_2.8.7        
##  [7] effects_4.2-0        carData_3.0-4        lattice_0.20-41     
## [10] future_1.21.0        rstan_2.21.2         StanHeaders_2.21.0-7
## [13] brms_2.14.4          Rcpp_1.0.5           simr_1.0.5          
## [16] lmerTest_3.1-3       lme4_1.1-26          Matrix_1.3-2        
## [19] knitr_1.30           egg_0.4.5            gridExtra_2.3       
## [22] forcats_0.5.0        stringr_1.4.0        dplyr_1.0.2         
## [25] purrr_0.3.4          readr_1.4.0          tidyr_1.1.2         
## [28] tibble_3.0.4         ggplot2_3.3.3        tidyverse_1.3.0     
## 
## loaded via a namespace (and not attached):
##   [1] tidyselect_1.1.0    htmlwidgets_1.5.3   grid_4.0.3         
##   [4] munsell_0.5.0       codetools_0.2-18    effectsize_0.4.1   
##   [7] statmod_1.4.35      DT_0.17             miniUI_0.1.1.1     
##  [10] withr_2.3.0         Brobdingnag_1.2-6   colorspace_2.0-0   
##  [13] rstudioapi_0.13     stats4_4.0.3        robustbase_0.93-7  
##  [16] bayesplot_1.8.0     listenv_0.8.0       emmeans_1.5.3      
##  [19] RLRsim_3.1-6        coda_0.19-4         parallelly_1.23.0  
##  [22] vctrs_0.3.6         generics_0.1.0      TH.data_1.0-10     
##  [25] xfun_0.20           R6_2.5.0            markdown_1.1       
##  [28] gamm4_0.2-6         projpred_2.0.2      assertthat_0.2.1   
##  [31] promises_1.1.1      scales_1.1.1        multcomp_1.4-15    
##  [34] nnet_7.3-14         gtable_0.3.0        globals_0.14.0     
##  [37] processx_3.4.5      sandwich_3.0-0      rlang_0.4.10       
##  [40] splines_4.0.3       broom_0.7.3         inline_0.3.17      
##  [43] yaml_2.2.1          reshape2_1.4.4      abind_1.4-5        
##  [46] modelr_0.1.8        threejs_0.3.3       crosstalk_1.1.1    
##  [49] backports_1.2.1     httpuv_1.5.5        rsconnect_0.8.16   
##  [52] tools_4.0.3         ellipsis_0.3.1      ggridges_0.5.3     
##  [55] plyr_1.8.6          base64enc_0.1-3     ps_1.5.0           
##  [58] prettyunits_1.1.1   zoo_1.8-8           haven_2.3.1        
##  [61] fs_1.5.0            survey_4.0          magrittr_2.0.1     
##  [64] data.table_1.13.6   openxlsx_4.2.3      colourpicker_1.1.0 
##  [67] reprex_0.3.0        mvtnorm_1.1-1       sjmisc_2.8.6       
##  [70] matrixStats_0.57.0  hms_1.0.0           shinyjs_2.0.0      
##  [73] mime_0.9            evaluate_0.14       xtable_1.8-4       
##  [76] shinystan_2.5.0     pbkrtest_0.5-0.1    rio_0.5.16         
##  [79] sjstats_0.18.1      readxl_1.3.1        ggeffects_1.0.1    
##  [82] rstantools_2.1.1    compiler_4.0.3      V8_3.4.0           
##  [85] crayon_1.3.4        minqa_1.2.4         htmltools_0.5.1    
##  [88] mgcv_1.8-33         later_1.1.0.1       RcppParallel_5.0.2 
##  [91] lubridate_1.7.9.2   DBI_1.1.0           sjlabelled_1.1.7   
##  [94] dbplyr_1.4.2        MASS_7.3-53         boot_1.3-25        
##  [97] cli_2.2.0           mitools_2.4         parallel_4.0.3     
## [100] insight_0.11.1      igraph_1.2.6        pkgconfig_2.0.3    
## [103] numDeriv_2016.8-1.1 foreign_0.8-81      binom_1.1-1        
## [106] xml2_1.3.2          dygraphs_1.1.1.6    estimability_1.3   
## [109] rvest_0.3.6         callr_3.5.1         digest_0.6.27      
## [112] parameters_0.10.1   fastGHQuad_1.0      rmarkdown_2.6      
## [115] cellranger_1.1.0    curl_4.3            shiny_1.5.0        
## [118] gtools_3.8.2        nloptr_1.2.2.2      lifecycle_0.2.0    
## [121] nlme_3.1-151        jsonlite_1.7.2      fansi_0.4.1        
## [124] pillar_1.4.7        loo_2.4.1           fastmap_1.0.1      
## [127] httr_1.4.2          plotrix_3.7-8       DEoptimR_1.0-8     
## [130] pkgbuild_1.2.0      survival_3.2-7      glue_1.4.2         
## [133] xts_0.12.1          bayestestR_0.8.0    zip_2.1.1          
## [136] shinythemes_1.1.2   iterators_1.0.13    stringi_1.5.3      
## [139] performance_0.6.1
# Read data ======================================
col_types <- cols(
  labid = col_factor(),
  subid = col_factor(),
  subid_unique = col_factor(),
  CDI.form = col_factor(),
  CDI.nwords = col_integer(),
  CDI.prop = col_number(),
  CDI.agerange = col_factor(),
  CDI.agedays = col_integer(),
  CDI.agemin = col_integer(),
  CDI.agemax = col_integer(),
  vocab_nwords = col_integer(),
  standardized.score.CDI = col_character(),
  standardized.score.CDI.num = col_number(),
  IDS_pref = col_number(),
  language = col_factor(),
  language_zone = col_factor(),
  CDI.error = col_logical(),
  Notes = col_character(),
  trial_order = col_factor(),
  method = col_factor(),
  age_days = col_integer(),
  age_mo = col_number(),
  age_group = col_factor(),
  nae = col_logical(),
  gender = col_factor(),
  second_session = col_logical()
)
data.total <- read_csv("data/02b_processed.csv", col_types = col_types)
# TODO: add saved results

Before moving on with the analysis, we have to ready the data by (a) checking for colinearity between z_age_months and CDI.z_age_months and correcting this if necessary, and (b) setting up the contrasts described in our data analysis.

1.1 Colinearity check

First, we run a Kappa test on the possibility of colinearity between z_age_months and CDI.z_age_months.

# Run kappa test
k.age_months <- model.matrix(~ z_age_months + CDI.z_age_months, data = data.total) %>%
  kappa(exact = T)

With a value of 3.1980728, we do not have a colinearity issue and can proceed with the data analysis as planned (The criteria of indicating colinearity is that kappa > 10).

1.2 Contrast Setups

We need gender as an effect-coded factor, and method as a deviation-coded factor. This is achieved in R by using the contr.sum() function with the number of levels for each factor. Notably, when subsetting the UK sample, only two levels of method out of the three in total were left.

# Set contrasts on the total dataset =============
contrasts(data.total$gender) <- contr.sum(2)
contrasts(data.total$method) <- contr.sum(3)

# Create sub-datasets, with contrasts ============
## NAE
data.nae <- data.total %>% subset(language_zone == "NAE") %>% droplevels()
contrasts(data.nae$gender) <- contr.sum(2)
contrasts(data.nae$method) <- contr.sum(3)

## UK
data.uk <- data.total %>% subset(language_zone == "British") %>% droplevels()
contrasts(data.uk$gender) <- contr.sum(2)
contrasts(data.uk$method) <- contr.sum(2) #note that UK sample has only 2 levels, so sum of zero contrasts set to 2 levels

## Other
data.other <- data.total %>% subset(language_zone == "Other") %>% droplevels()
contrasts(data.other$gender) <- contr.sum(2)
contrasts(data.other$method) <- contr.sum(3)

2 Descriptive Statistics

We first assess the amount of data we have overall per condition and their shape overall.

data.total %>%
  group_by(language_zone, CDI.agerange, method, gender) %>%
  summarise(N = n(), age = mean(CDI.agedays), sd = sd(CDI.agedays)) %>%
  kable()
## `summarise()` regrouping output by 'language_zone', 'CDI.agerange', 'method' (override with `.groups` argument)
language_zone CDI.agerange method gender N age sd
British 18 singlescreen F 24 551.2083 10.741950
British 18 singlescreen M 24 550.5000 13.497182
British 18 eyetracking F 9 548.5556 9.139353
British 18 eyetracking M 11 547.8182 10.147100
British 24 singlescreen F 18 741.6111 13.128948
British 24 singlescreen M 15 745.0667 15.549307
British 24 eyetracking F 13 731.5385 12.162890
British 24 eyetracking M 5 737.8000 8.228001
Other 18 singlescreen F 11 541.5455 7.160498
Other 18 singlescreen M 13 544.3077 6.663140
Other 18 eyetracking F 26 556.0769 14.133430
Other 18 eyetracking M 30 560.8333 16.128970
Other 18 hpp F 38 549.0526 10.812774
Other 18 hpp M 38 555.1579 13.664961
Other 24 singlescreen F 10 731.3000 14.802777
Other 24 singlescreen M 12 721.1667 13.888081
Other 24 eyetracking F 28 730.1786 9.996494
Other 24 eyetracking M 25 730.1600 7.711896
Other 24 hpp F 45 730.7333 17.357733
Other 24 hpp M 35 730.5143 15.849237
NAE 18 singlescreen F 19 554.9474 20.780530
NAE 18 singlescreen M 14 581.2143 24.925052
NAE 18 eyetracking F 1 549.0000 NA
NAE 18 hpp F 64 557.9375 22.801333
NAE 18 hpp M 66 556.1667 22.762035
NAE 24 singlescreen F 23 737.1739 26.604258
NAE 24 singlescreen M 20 739.6000 21.352678
NAE 24 eyetracking F 2 756.5000 34.648232
NAE 24 eyetracking M 1 745.0000 NA
NAE 24 hpp F 58 731.6897 28.149480
NAE 24 hpp M 65 726.2769 27.133068

More detailed information about Descriptive Statistics

#number of lab
data.total %>% 
  select(labid, language_zone) %>% 
  unique() %>% 
  group_by(language_zone) %>% 
  count()
## # A tibble: 3 x 2
## # Groups:   language_zone [3]
##   language_zone     n
##   <fct>         <int>
## 1 British           4
## 2 Other             8
## 3 NAE               9
data.total %>% 
  group_by(language_zone, CDI.agerange) %>% 
  summarize(N = n())
## `summarise()` regrouping output by 'language_zone' (override with `.groups` argument)
## # A tibble: 6 x 3
## # Groups:   language_zone [3]
##   language_zone CDI.agerange     N
##   <fct>         <fct>        <int>
## 1 British       18              68
## 2 British       24              51
## 3 Other         18             156
## 4 Other         24             155
## 5 NAE           18             164
## 6 NAE           24             169
# age range in each age group and language zone
data.total %>% 
  select(subid, language_zone, CDI.agedays, CDI.agerange) %>% 
  unique() %>% 
  group_by(language_zone, CDI.agerange) %>% 
  summarize(age_min = (min(CDI.agedays)/365.25*12),
            age_max = (max(CDI.agedays)/365.25*12))
## `summarise()` regrouping output by 'language_zone' (override with `.groups` argument)
## # A tibble: 6 x 4
## # Groups:   language_zone [3]
##   language_zone CDI.agerange age_min age_max
##   <fct>         <fct>          <dbl>   <dbl>
## 1 British       18              17.4    19.2
## 2 British       24              23.4    25.1
## 3 Other         18              17.4    19.5
## 4 Other         24              22.5    25.2
## 5 NAE           18              16.1    20.0
## 6 NAE           24              22.1    25.9

We then assess the data per lab in terms of sample size and CDI score (vocabulary size, for consistency between language zones).

by_lab <- data.total %>%
  group_by(labid, language_zone, CDI.agerange) %>%
  mutate(tested = n_distinct(subid_unique)) %>%
  select(labid, language_zone, CDI.agerange, tested, vocab_nwords) %>%
  nest(scores = vocab_nwords) %>%
  mutate(model = map(scores, ~ lm(vocab_nwords ~ 1, data = .x)),
         ci = map(model, confint)) %>%
  transmute(tested = tested,
            mean = map_dbl(model, ~ coefficients(.x)[[1]]),
            ci_lower = map_dbl(ci, 1),
            ci_upper = map_dbl(ci, 2)) %>%
  arrange(language_zone) %>%
  rownames_to_column()

# TODO: find a way to group by language zone?
ggplot(by_lab,
       aes(x = labid, colour = language_zone,
           y = mean, ymin = ci_lower, ymax = ci_upper)) + 
  geom_linerange() + 
  geom_point(aes(size = tested)) + 
  facet_grid(cols = vars(CDI.agerange), scales = "free") + coord_flip(ylim = c(0, 500)) +
  xlab("Lab") + ylab("Vocabulary size") +
  scale_colour_brewer(palette = "Dark2", name = "Language zone",
                      breaks = c("British", "NAE", "Other")) +
  scale_size_continuous(name = "N", breaks = function(x) c(min(x), mean(x), max(x))) +
  theme(legend.position = "bottom")

3 Sample Theory Based Statistics

3.1 Simple Correlation

First, we want to assess quickly if there is a direct correlation between IDS preference and CDI score, computing a Pearson#’s product-moment correlation. We use standardized CDI scores for the North American sample, and raw scores for the British sample. [NOTE FROM MELANIE: I’m not sure this makes sense to do with the raw UK scores - since we’re collapsing across 18 and 24 month data and the data aren#’t standardized by age.].

# Statistics =====================================
## North American Sample
test.pearson.nae <- cor.test(data.nae$IDS_pref,
                             data.nae$z_standardized_CDI,
                             alternative = "two.sided", method = "pearson")

test.pearson.nae
## 
##  Pearson's product-moment correlation
## 
## data:  data.nae$IDS_pref and data.nae$z_standardized_CDI
## t = -0.91847, df = 305, p-value = 0.3591
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.16349833  0.05977293
## sample estimates:
##         cor 
## -0.05251901
## UK Sample
test.pearson.uk <- cor.test(data.uk$IDS_pref,
                            data.uk$z_vocab_nwords,
                            alternative = "two.sided", method = "pearson")

test.pearson.uk
## 
##  Pearson's product-moment correlation
## 
## data:  data.uk$IDS_pref and data.uk$z_vocab_nwords
## t = 1.0327, df = 117, p-value = 0.3039
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.08643398  0.27040987
## sample estimates:
##        cor 
## 0.09504018

Plots for correlation

## North American Sample
### Get correlation value for annotation
cor_text <- "paste(italic(R)^2, \" =\")"
cor_value <- round(test.pearson.nae$estimate, 3)

### Build plot
plot.pearson.nae <- data.nae %>%
  ggplot(aes(x = IDS_pref,
             y = standardized.score.CDI.num)) +
  xlab("IDS preference") + ylab("Standardized CDI score") +
  geom_point() +
  geom_smooth(method = lm) +
  annotate("text", x = -.9, y = 51, parse = T, size = 4,
           label = paste(cor_text, cor_value, sep = "~"))

## UK Sample
cor_value <- round(test.pearson.uk$estimate, 3)
plot.pearson.uk <- data.uk %>%
  ggplot(aes(x = IDS_pref,
             y = vocab_nwords)) +
  xlab("IDS preference") + ylab("Vocabulary size (in words)") +
  geom_point() +
  geom_smooth(method = lm) +
  annotate("text", x = .8, y = 150, parse = T, size = 4,
           label = paste(cor_text, cor_value, sep = "~"))

# Global plot
plot.pearson <- ggarrange(plot.pearson.nae, plot.pearson.uk, ncol = 2)
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 26 rows containing non-finite values (stat_smooth).
## Warning: Removed 26 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'

plot.pearson

ggsave("plots/corr_plot.pdf", plot.pearson,
       units = "mm", width = 180, height = 100, dpi = 1000)

We see no obvious direct link between IDS prefernce and CDI score here. However, an effect might appear once we take into account various factors that might interact with IDS preference and/or CDI score. We can also first enhance these plots with information about the age group at which infants were tested (18- or 24-month-old), using vocabulary size to better compare the NAE and UK samples.

plot.age_group <- data.total %>%
  subset(language_zone != "Other") %>%
  droplevels() %>%
  ggplot(aes(x = IDS_pref,
             y = vocab_nwords,
             colour = CDI.agerange)) +
  facet_wrap(vars(language_zone),
             labeller = as_labeller(c("British" = "UK samples",
                                      "NAE" = "North Amercian samples"))) +
  xlab("Standardized IDS prefence score") + ylab("Vocabulary size (in words)") + theme(legend.position = "top") +
  geom_point() +
  geom_smooth(method = lm) +
  scale_colour_brewer(palette = "Dark2", name = "Age group",
                      breaks = c("18", "24"),
                      labels = c("18mo", "24m"))
ggsave("plots/scatter_age.pdf", plot.age_group,
       units = "mm", width = 180, height = 100, dpi = 1000)
## `geom_smooth()` using formula 'y ~ x'
(plot.age_group)
## `geom_smooth()` using formula 'y ~ x'

3.2 Mixed-Effects Model by Language Zone

Here, we run a mixed-effects model including only theoretically motivated effects, as described in the pre-registration. We start with the full model bellow, simplifying the random effects structure until it converges.

3.2.1 NAE full model

# Run models =====================================
## NAE

data.nae$centered_IDS_pref <- scale(data.nae$IDS_pref, center = T, scale = F)

lmer.full.nae <- lmer(standardized.score.CDI.num ~ CDI.z_age_months + gender +
                        z_age_months + method + centered_IDS_pref +
                        centered_IDS_pref:method + centered_IDS_pref:CDI.z_age_months + centered_IDS_pref:z_age_months +
                        (1 | labid) + (1 | subid_unique),
                      data = data.nae)

summary(lmer.full.nae)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months +  
##     method + centered_IDS_pref + centered_IDS_pref:method + centered_IDS_pref:CDI.z_age_months +  
##     centered_IDS_pref:z_age_months + (1 | labid) + (1 | subid_unique)
##    Data: data.nae
## 
## REML criterion at convergence: 2806.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2240 -0.4550 -0.0619  0.4779  2.5238 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  subid_unique (Intercept) 531.62   23.06   
##  labid        (Intercept)  37.34    6.11   
##  Residual                 246.81   15.71   
## Number of obs: 307, groups:  subid_unique, 211; labid, 8
## 
## Fixed effects:
##                                     Estimate Std. Error        df t value
## (Intercept)                         32.47880    6.47774  49.79324   5.014
## CDI.z_age_months                     0.04723    0.35956 121.44454   0.131
## gender1                             -1.51119    1.88476 193.54899  -0.802
## z_age_months                        -0.01368    0.62999  69.41345  -0.022
## method1                              6.62448    6.56943 122.84082   1.008
## method2                            -13.97899   11.88500 177.54529  -1.176
## centered_IDS_pref                  -34.17431   24.81975 210.14326  -1.377
## method1:centered_IDS_pref           28.35879   25.46036 210.29032   1.114
## method2:centered_IDS_pref          -59.58231   49.55101 209.34368  -1.202
## CDI.z_age_months:centered_IDS_pref   0.41746    1.10252 123.90420   0.379
## z_age_months:centered_IDS_pref       2.01470    2.07816 189.23346   0.969
##                                    Pr(>|t|)    
## (Intercept)                        7.14e-06 ***
## CDI.z_age_months                      0.896    
## gender1                               0.424    
## z_age_months                          0.983    
## method1                               0.315    
## method2                               0.241    
## centered_IDS_pref                     0.170    
## method1:centered_IDS_pref             0.267    
## method2:centered_IDS_pref             0.231    
## CDI.z_age_months:centered_IDS_pref    0.706    
## z_age_months:centered_IDS_pref        0.334    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) CDI.z__ gendr1 z_g_mn methd1 methd2 c_IDS_ m1:_ID m2:_ID
## CDI.z_g_mnt -0.058                                                         
## gender1     -0.031  0.022                                                  
## z_age_mnths  0.041 -0.037  -0.021                                          
## method1     -0.716  0.000   0.021  0.105                                   
## method2      0.880 -0.030  -0.035  0.017 -0.877                            
## cntrd_IDS_p  0.351  0.006   0.038 -0.040 -0.342  0.380                     
## mthd1:_IDS_ -0.330 -0.028  -0.027  0.017  0.356 -0.379 -0.927              
## mthd2:_IDS_  0.344  0.021   0.060 -0.008 -0.356  0.386  0.966 -0.971       
## CDI.__:_IDS -0.011  0.012  -0.019 -0.019 -0.029  0.008 -0.049  0.027 -0.047
## z_g_m:_IDS_ -0.041  0.026   0.103  0.047 -0.037  0.022  0.028 -0.086  0.132
##             CDI.__:
## CDI.z_g_mnt        
## gender1            
## z_age_mnths        
## method1            
## method2            
## cntrd_IDS_p        
## mthd1:_IDS_        
## mthd2:_IDS_        
## CDI.__:_IDS        
## z_g_m:_IDS_ -0.071
# robust_lmer.full.nae <- robustlmm::rlmer(standardized.score.CDI.num ~ CDI.z_age_months + gender +
#                         z_age_months + method + centered_IDS_pref +
#                         centered_IDS_pref:method + centered_IDS_pref:CDI.z_age_months + centered_IDS_pref:z_age_months +
#                         (1 | labid),
#                       data = data.nae)
# 
# 
# summary(robust_lmer.full.nae) #this model is used to see if we can meet some statistical assumption better but we decided to use the original model as the inferential statistical results are consistent

full.nae_pvalue <- anova(lmer.full.nae) %>% 
  as_tibble(rownames = "Parameter") #this gives us the Type III p values

# ==========
# Sequentially removed random effects:
# IDS_pref:z_age_months
# IDS_pref:CDI.z_age_months
# IDS_pref:method
# IDS_pref
# method
# z_age_months
# gender
# CDI.z_age_months
#==========

3.2.1.1 (Optional) Checking mixed-model assumptions. We will check the following:

  1. Linearlity
  2. Normality of the residuals
  3. Homoscedasticity of residuals
  4. No autocorrelation
  5. No multicollinearity
#First, check linearity
# data.nae$resid <- residuals(lmer.full.nae)
# 
# plot(data.nae$resid, data.nae$standardized.score.CDI)

#Second, check normality
plot_model(lmer.full.nae, type = 'diag') ## we do have right-skewed normality of residuals
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.18
## Current Matrix version is 1.3.2
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
## [[1]]
## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]
## [[2]]$subid_unique
## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]$labid
## `geom_smooth()` using formula 'y ~ x'

## 
## 
## [[3]]

## 
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

#Third, check autocorrelation
re_run_lme.full.nae <- nlme::lme(standardized.score.CDI.num ~ CDI.z_age_months + gender +
                        z_age_months + method + centered_IDS_pref +
                        centered_IDS_pref:method + centered_IDS_pref:CDI.z_age_months + centered_IDS_pref:z_age_months, random = ~1 | labid,
                        method = "REML",
                      data = data.nae, na.action = na.exclude)

plot(nlme::ACF(re_run_lme.full.nae, resType = "normalized")) #there is no sign for autocorrelation

#Lastly, check multi-collinearity
car::vif(lmer.full.nae) #we do see a multicollineartiy for the IDS preference variable, even though we have centered the IDS preference score. It is probably related to the number the participating labs (as this is the group level that we are controlling) and how we entered interaction between IDS preference and other variables (that lack variability in the current sample). We need to keep IDS preference in the model as exploring the relationship between IDS preference and CDI score is the key research question in the paper.
##                                         GVIF Df GVIF^(1/(2*Df))
## CDI.z_age_months                    1.009314  1        1.004646
## gender                              1.039338  1        1.019479
## z_age_months                        1.096811  1        1.047288
## method                              1.284501  2        1.064593
## centered_IDS_pref                  19.298596  1        4.393017
## method:centered_IDS_pref           20.886393  2        2.137794
## CDI.z_age_months:centered_IDS_pref  1.014832  1        1.007389
## z_age_months:centered_IDS_pref      1.292422  1        1.136847

3.2.2 UK full model

lmer.full.uk <- lmer(vocab_nwords ~ CDI.z_age_months + gender +
                       z_age_months + method + IDS_pref +
                       IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
                       (1 | labid) + (1 | subid_unique),
                     data = data.uk) 

summary(lmer.full.uk)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method +  
##     IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months +  
##     IDS_pref:z_age_months + (1 | labid) + (1 | subid_unique)
##    Data: data.uk
## 
## REML criterion at convergence: 1326
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.78175 -0.38449  0.01193  0.37856  1.94702 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  subid_unique (Intercept) 5412.10  73.567  
##  labid        (Intercept)   33.24   5.765  
##  Residual                 2637.21  51.354  
## Number of obs: 119, groups:  subid_unique, 88; labid, 4
## 
## Fixed effects:
##                           Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)                174.753     11.544   3.754  15.138 0.000167 ***
## CDI.z_age_months            28.905      1.964  44.974  14.714  < 2e-16 ***
## gender1                     19.043      9.820  78.517   1.939 0.056074 .  
## z_age_months                -5.065      3.455  70.444  -1.466 0.147086    
## method1                    -11.249     13.255   5.502  -0.849 0.431430    
## IDS_pref                     9.818     26.555  82.944   0.370 0.712523    
## method1:IDS_pref            -7.789     29.223  87.125  -0.267 0.790454    
## CDI.z_age_months:IDS_pref    1.573      5.726  50.813   0.275 0.784675    
## z_age_months:IDS_pref        1.154      8.805  89.735   0.131 0.896024    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) CDI.z__ gendr1 z_g_mn methd1 IDS_pr m1:IDS CDI.__:
## CDI.z_g_mnt  0.120                                                   
## gender1      0.044 -0.069                                            
## z_age_mnths -0.011 -0.099  -0.019                                    
## method1     -0.395 -0.074  -0.090  0.490                             
## IDS_pref    -0.360 -0.074  -0.112  0.033  0.199                      
## mthd1:IDS_p  0.188  0.094   0.007 -0.117 -0.303 -0.203               
## CDI.__:IDS_ -0.099 -0.328   0.056  0.064  0.094  0.090 -0.136        
## z_g_mn:IDS_ -0.025  0.079  -0.244 -0.374 -0.099 -0.112  0.431 -0.228
full.uk_pvalue <- anova(lmer.full.uk) %>% 
  as_tibble(rownames = "Parameter") #this gives us the Type III p values

#==========
# Sequentially removed random effects:
# IDS_pref:z_age_months
# IDS_pref:CDI.z_age_months
# IDS_pref:method
# IDS_pref
# method
# z_age_months
# gender
# CDI.z_age_months

3.2.2.1 (Optional) Checking mixed-model assumptions. We will check the following:

  1. Linearlity
  2. Normality of the residuals
  3. Homoscedasticity of residuals
  4. No autocorrelation
  5. No multicollinearity
#First, check linearity. The plot looks linear
data.uk$resid <- residuals(lmer.full.uk)
 
plot(data.uk$resid, data.uk$vocab_nwords)

#Second, check normality
plot_model(lmer.full.uk, type = 'diag') ## we do have right-skewed normality of residuals
## [[1]]
## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]
## [[2]]$subid_unique
## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]$labid
## `geom_smooth()` using formula 'y ~ x'

## 
## 
## [[3]]

## 
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

#Third, check autocorrelation
re_run_lme.full.uk <- nlme::lme(vocab_nwords ~ CDI.z_age_months + gender +
                       z_age_months + method + IDS_pref +
                       IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months, random = ~1 | labid,
                        method = "REML",
                      data = data.nae, na.action = na.exclude)

plot(nlme::ACF(re_run_lme.full.uk, resType = "normalized")) #there is no sign for autocorrelation

#Lastly, check multi-collinearity
car::vif(lmer.full.uk) #no problem for multicollinearlity
##          CDI.z_age_months                    gender              z_age_months 
##                  1.142786                  1.125370                  1.678292 
##                    method                  IDS_pref           method:IDS_pref 
##                  1.592427                  1.099496                  1.457526 
## CDI.z_age_months:IDS_pref     z_age_months:IDS_pref 
##                  1.188763                  1.737367

We now want to check the statistical power of significant effects, and discard any models with significant effects that do not reach 80% power. This however leads to too many warnings of singularity issues on the model updates inherent to the simr power simulations, hence we cannot obtain satisfactory power estimates as pre-registered.

AST: Note that we don’t have any IV(s) that turned out to be significant in the Full NAE model. So we won’t run the power analysis check. For the UK full model, there are two statistically significant IV: CDI_age and gender. The post hoc power check suggested that we have high power in detecting the effect of CDI_age but not gender. Note that gender has a smaller effect size to begin with, so this may partially explain why we have less power in detecting it in the model. As there can be a number of different factors that determines the posthoc power, we decided not to remove gender in the model based on posthoc power analysis check.

check_pwr_uk_cdi_age <- simr::powerSim(lmer.full.uk, test = fixed("CDI.z_age_months", method = "z"), seed = 2, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into

check_pwr_uk_cdi_age

check_pwr_uk_gender <- simr::powerSim(lmer.full.uk, test = fixed("gender", method = "z"), seed = 3, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into

check_pwr_uk_gender

3.2.3 Combined Sample

For this combined analysis, we first need to restrain the age range for the NAE sample (previously ±2 months, now ±1.5 months).

# Create dataset with British and NAE only
data.uk_nae <- data.total %>%
  subset(language_zone %in% c("British", "NAE")) %>%
  mutate(CDI.agemin = ifelse(language_zone == "NAE",
                             CDI.agemin + round(.5*365.25/12),
                             CDI.agemin),
         CDI.agemax = ifelse(language_zone == "NAE",
                             CDI.agemax - round(.5*365.25/12),
                             CDI.agemax)) %>%
  subset(!(CDI.agedays < CDI.agemin | CDI.agedays > CDI.agemax)) %>%
  droplevels()
# Create contrasts for analysis
contrasts(data.uk_nae$gender) <- contr.sum(2)
contrasts(data.uk_nae$method) <- contr.sum(3)
contrasts(data.uk_nae$language_zone) <- contr.sum(2)

We can then run the planned combined analysis adding the main effect and interactions of language_zone.

lmer.full.uk_nae <- lmer(CDI.prop ~ CDI.z_age_months + language_zone + gender +
                           z_age_months + method + IDS_pref + IDS_pref:language_zone +
                           IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
                           (1 + CDI.z_age_months | labid) + (1 | subid_unique),
                         data = data.uk_nae)

summary(lmer.full.uk_nae)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CDI.prop ~ CDI.z_age_months + language_zone + gender + z_age_months +  
##     method + IDS_pref + IDS_pref:language_zone + IDS_pref:method +  
##     IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 +  
##     CDI.z_age_months | labid) + (1 | subid_unique)
##    Data: data.uk_nae
## 
## REML criterion at convergence: -61.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.83614 -0.46389 -0.04311  0.41335  2.60155 
## 
## Random effects:
##  Groups       Name             Variance  Std.Dev. Corr
##  subid_unique (Intercept)      0.0267932 0.16369      
##  labid        (Intercept)      0.0003603 0.01898      
##               CDI.z_age_months 0.0009806 0.03131  0.44
##  Residual                      0.0184535 0.13584      
## Number of obs: 424, groups:  subid_unique, 302; labid, 13
## 
## Fixed effects:
##                             Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                 0.320830   0.017878  12.110147  17.945 4.31e-10 ***
## CDI.z_age_months            0.043700   0.009272  13.321953   4.713  0.00038 ***
## language_zone1              0.087738   0.021795  10.500676   4.026  0.00219 ** 
## gender1                     0.032723   0.011939 270.491415   2.741  0.00654 ** 
## z_age_months               -0.002774   0.004295 164.183977  -0.646  0.51919    
## method1                    -0.007741   0.024041  21.324230  -0.322  0.75061    
## method2                     0.009802   0.035669  15.572064   0.275  0.78708    
## IDS_pref                    0.032520   0.041549 295.409651   0.783  0.43445    
## language_zone1:IDS_pref     0.036505   0.056770 318.273301   0.643  0.52066    
## method1:IDS_pref           -0.040507   0.053921 317.820649  -0.751  0.45307    
## method2:IDS_pref           -0.018199   0.087818 306.828123  -0.207  0.83596    
## CDI.z_age_months:IDS_pref   0.009728   0.008057 172.834552   1.207  0.22889    
## z_age_months:IDS_pref      -0.001722   0.011708 278.006432  -0.147  0.88318    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
combined.full.uk_nae_pvalue <- anova(lmer.full.uk_nae) %>% 
  as_tibble(rownames = "Parameter") #this gives us the Type III p values

#==========
# Sequentially removed random effects:
# IDS_pref:z_age_months
# IDS_pref:CDI.z_age_months
# IDS_pref:method
# IDS_pref:language_zone
# IDS_pref
# method
# z_age_months
# gender
# language_zone
#==========

3.2.3.1 (Optional) Checking mixed-model assumptions

  1. Linearlity
  2. Normality of the residuals
  3. Homoscedasticity of residuals
  4. No autocorrelation
  5. No multicollinearity
#First, check linearity. The plot looks linear
data.uk_nae$resid <- residuals(lmer.full.uk_nae)
 
plot(data.uk_nae$resid, data.uk_nae$CDI.prop)

#Second, check normality
plot_model(lmer.full.uk_nae, type = 'diag') ## we do have right-skewed normality of residuals
## [[1]]
## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]
## [[2]]$subid_unique
## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]$labid
## `geom_smooth()` using formula 'y ~ x'

## 
## 
## [[3]]

## 
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

#Third, check autocorrelation
re_run_lme.full.uk_nae <- nlme::lme(CDI.prop ~ CDI.z_age_months + language_zone + gender +
                           z_age_months + method + IDS_pref + IDS_pref:language_zone +
                           IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months,
                           random = ~CDI.z_age_months  | labid,
                        method = "REML",
                      data = data.uk_nae, na.action = na.exclude)

plot(nlme::ACF(re_run_lme.full.uk_nae, resType = "normalized")) #there is no sign for autocorrelation

#Lastly, check multi-collinearity
car::vif(lmer.full.uk_nae) #no problem for multicollinearlity
##                               GVIF Df GVIF^(1/(2*Df))
## CDI.z_age_months          1.019950  1        1.009926
## language_zone             2.325115  1        1.524833
## gender                    1.036435  1        1.018055
## z_age_months              1.474609  1        1.214335
## method                    3.284643  2        1.346239
## IDS_pref                  1.447535  1        1.203135
## language_zone:IDS_pref    2.971184  1        1.723712
## method:IDS_pref           3.985056  2        1.412891
## CDI.z_age_months:IDS_pref 1.053846  1        1.026570
## z_age_months:IDS_pref     1.430296  1        1.195950

We then compute \(p\)-values, but leave out power estimates for those \(p\)-values as above. Again, we have a lot of singular fit issues for the power checks and decided not to remove parameters based on posthoc power analysis.

check_pwr_combined_cdi_age <- simr::powerSim(lmer.full.uk_nae, test = fixed("CDI.z_age_months", method = "z"), seed = 2, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into

check_pwr_combined_cdi_age

check_pwr_combined_lang_zone <- simr::powerSim(lmer.full.uk_nae, test = fixed("language_zone", method = "z"), seed = 3, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into

check_pwr_combined_lang_zone

check_pwr_combined_gender <- simr::powerSim(lmer.full.uk_nae, test = fixed("gender", method = "z"), seed = 3, nsim = 1000, alpha = 0.05) #specify that Gender is the ixed effect that we are looking into

check_pwr_combined_gender

4 Bayesian Statistics

Given that the `frequentist#’ tests above did not discover any expected significant effects, we now run a Bayesian analysis, specifically Bayes factors (or \(b\)-values) from model comparisons, to determine whether our data provide evidence for a true null effect or merely fail to provide evidence either way.

We first run models on the separate UK and NAE samples.

First of all, we need to set up the priors for the Bayeisan analysis. In the RR, we said that we would use a Cauchy distribution for priors over fixed effects with normally distributed priors over random effect parameters."

4.1 Set up priors

Set up a really weak and not too informative prior

prior <-c(set_prior("normal(10, 1)", class = "Intercept"), #we can't have negative CDI
              set_prior("cauchy(0, 1)", class = "b"), #all IVs have the same prior..
              set_prior("normal(0, 1)", class = "sd")) 

prior_combined_mod <-c(set_prior("normal(10, 1)", class = "Intercept"),
              set_prior("cauchy(0, 1)", class = "b"), #all IVs have the same prior..
              set_prior("normal(0, 1)", class = "sd"),
              set_prior("lkj(2)", class = "L")) #adding this as there is a random slope and we need a correlation matrix prior between slope and intercept

Set up another prior that mainly changed the meaan for the fixed effect and see if there is big difference? Note: numbers changed but not to the extent that change BF conclusion.

prior_big <-c(set_prior("normal(10, 1)", class = "Intercept"), #we can't have negative CDI
              set_prior("cauchy(10, 5)", class = "b"), #all IVs have the same prior..
              set_prior("normal(0, 1)", class = "sd")) 

prior_combined_big <-c(set_prior("normal(10, 1)", class = "Intercept"),
              set_prior("cauchy(10, 5)", class = "b"), #all IVs have the same prior..
              set_prior("normal(0, 1)", class = "sd"),
              set_prior("lkj(2)", class = "L")) #adding this as there is a random slope and we need a correlation matrix prior between slope and intercept

We also try more informative prior.

4.2 NAE info prior

prior_nae_1 <-c(set_prior("normal(50, 10)", class = "Intercept"), #note that nae's DV is percentile score, but I will assume at 18 months, the average is at 50 percentile
              set_prior("cauchy(5, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 5 % percentile increase?
              set_prior("cauchy(-15, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males
              set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
               set_prior("normal(0, 1)", class = "sd")) 

prior_nae_null_1 <-c(set_prior("normal(50, 10)", class = "Intercept"), #note that nae's DV is percentile score, but I will assume at 18 months, the average is at 50 percentile
              set_prior("cauchy(5, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 5 % percentile increase?
              set_prior("cauchy(-15, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
               set_prior("normal(0, 1)", class = "sd")) 

prior_nae_2 <-c(set_prior("normal(50, 10)", class = "Intercept"), #note that nae's DV is percentile score, but I will assume at 18 months, the average is at 50 percentile
              set_prior("cauchy(5, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 5 % percentile increase?
              set_prior("cauchy(-15, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males
              set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
              #note all interaction priors are weak as we don't really have much info to form a prior
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
              set_prior("cauchy(0, 1)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method1:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method2:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months:IDS_pref"),
              set_prior("normal(0, 1)", class = "sd")) 

prior_nae_null_2 <-c(set_prior("normal(50, 10)", class = "Intercept"), #note that nae's DV is percentile score, but I will assume at 18 months, the average is at 50 percentile
              set_prior("cauchy(5, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 5 % percentile increase?
              set_prior("cauchy(-15, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males
              set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores
               #note all interaction priors are weak as we don't really have much info to form a prior
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
              set_prior("cauchy(0, 1)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method1:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method2:IDS_pref"),
              set_prior("normal(0, 1)", class = "sd")) 

4.3 UK info prior

prior_uk_1 <-c(set_prior("normal(40, 20)", class = "Intercept"), #note that UK's DV is raw score, based on wordbank (oxford CDI), the average (50 percentile) is around 41.6 
              set_prior("cauchy(15, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase (proxy from wordbank)
              set_prior("cauchy(-10, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males
              set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
              set_prior("normal(0, 1)", class = "sd")) 

prior_uk_null_1 <-c(set_prior("normal(40, 20)", class = "Intercept"), #note that UK's DV is raw score, based on wordbank (oxford CDI), the average (50 percentile) is around 41.6 
              set_prior("cauchy(15, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase (proxy from wordbank)
              set_prior("cauchy(-10, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
              set_prior("normal(0, 1)", class = "sd")) 

prior_uk_2 <-c(set_prior("normal(40, 20)", class = "Intercept"), #note that UK's DV is raw score, based on wordbank (oxford CDI), the average (50 percentile) is around 41.6 
              set_prior("cauchy(15, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase (proxy from wordbank)
              set_prior("cauchy(-10, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males
              set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
               #note all interaction priors are weak as we don't really have much info to form a prior
              set_prior("cauchy(0, 1)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method1:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months:IDS_pref"),
              set_prior("normal(0, 1)", class = "sd")) 

prior_uk_null_2 <-c(set_prior("normal(40, 20)", class = "Intercept"), #note that UK's DV is raw score, based on wordbank (oxford CDI), the average (50 percentile) is around 41.6 
              set_prior("cauchy(15, 1)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase (proxy from wordbank)
              set_prior("cauchy(-10, 5)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males
              set_prior("cauchy(0, 1)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 1)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 1)", class = "b", coef = "z_age_months"),
               #note all interaction priors are weak as we don't really have much info to form a prior
              set_prior("cauchy(0, 1)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
              set_prior("cauchy(0, 1)", class = "b", coef = "method1:IDS_pref"),
              set_prior("normal(0, 1)", class = "sd")) 

4.4 Combined info prior

# Note that NAE CDI WS total score is 680 and UK CDI total score is 416. 

# NAE average scores for 50th percentile is 76.3, so it is 76.3/680 ~ 11% and UK CDI average scores for 50th percentile is 41.6, so it is 41.6/416 ~ 10%

prior_combined_full_info <-c(set_prior("normal(0.1, 0.05)", class = "Intercept"), # given the above info, we expect on average, proportion score is 10% 
              set_prior("cauchy(0.03, 0.01)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase in UK cdi and lead to 28 words increase in NAE (proxy from wordbank). In general, 15/416 = 3% and 28/680 = 4% 
              set_prior("cauchy(-0.03, 0.01)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males, 11/416 ~ 2.6% and Wordbank NAE suggested females knows approximately 27 words than males, 27/680 ~ 3.9%
               set_prior("cauchy(0, 0.01)", class = "b", coef = "language_zone1"), #weak prior as don't know the effect of language zone to CDI score 
              set_prior("cauchy(0, 0.01)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 0.01)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 0.01)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 0.01)", class = "b", coef = "z_age_months"),
               #note all interaction priors are weak as we don't really have much info to form a prior
              set_prior("cauchy(0, 0.01)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
              set_prior("cauchy(0, 0.01)", class = "b", coef = "language_zone1:IDS_pref"),
              set_prior("cauchy(0, 0.01)", class = "b", coef = "method1:IDS_pref"),
              set_prior("cauchy(0, 0.01)", class = "b", coef = "method2:IDS_pref"),
              set_prior("cauchy(0, 0.01)", class = "b", coef = "z_age_months:IDS_pref"),
              set_prior("normal(0, 0.01)", class = "sd"),
              set_prior("lkj(2)", class = "L")) #adding this as there is a random slope and we need a correlation matrix prior between slope and intercept

prior_combined_null_info <-c(set_prior("normal(0.1, 0.05)", class = "Intercept"), # given the above info, we expect on average, proportion score is 10% 
              set_prior("cauchy(0.03, 0.01)", class = "b", coef = "CDI.z_age_months"), #assume increase in 1 month may lead to 15 words increase in UK cdi and lead to 28 words increase in NAE (proxy from wordbank). In general, 15/416 = 3% and 28/680 = 4% 
              set_prior("cauchy(-0.03, 0.01)", class = "b", coef = "gender1"), # female is the baseline, normally females know more words than males. Wordbank oxford CDI suggested at 18 months (50 percentile), female knows approximately 11 words more than males, 11/416 ~ 2.6% and Wordbank NAE suggested females knows approximately 27 words than males, 27/680 ~ 3.9%
               set_prior("cauchy(0, 0.01)", class = "b", coef = "language_zone1"), #weak prior as don't know the effect of language zone to CDI score 
              set_prior("cauchy(0, 0.01)", class = "b", coef = "IDS_pref"), #weak prior as don't know the effect of IDS preference to CDI score
              set_prior("cauchy(0, 0.01)", class = "b", coef = "method1"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 0.01)", class = "b", coef = "method2"), #weak prior as methods related to the MB1 study shouldn't influence CDI scores 
              set_prior("cauchy(0, 0.01)", class = "b", coef = "z_age_months"),
               #note all interaction priors are weak as we don't really have much info to form a prior
              set_prior("cauchy(0, 0.01)", class = "b", coef = "CDI.z_age_months:IDS_pref"),
              set_prior("cauchy(0, 0.01)", class = "b", coef = "method1:IDS_pref"),
              set_prior("cauchy(0, 0.01)", class = "b", coef = "method2:IDS_pref"),
              set_prior("cauchy(0, 0.01)", class = "b", coef = "z_age_months:IDS_pref"),
              set_prior("normal(0, 0.01)", class = "sd"),
              set_prior("lkj(2)", class = "L")) #adding this as there is a random slope and we need a correlation matrix prior between slope and intercept

In the following analysis, we follow what we have proposed in the RR about calculating the Bayes factor for the three variables of interest: 1) The main effect IDS preference on CDI scores 2) The interaction effect between IDS preference and IDS test age on CDI scores 3) The interaction effect between language_zone (i.e., dialect in RR) and IDS preference on CDI scores

For the first two variables, we run it on the NAE and UK model separately. For the third variable, we run in on the combined UK & NAE model.

To calcualte Bayes factor for the main effect IDS preference on CDI scores, we need to two models: a full model that contains the main effect of variables of interest, and the other null model that does not contain the main effect of interest. As all the interaction terms in the models contains IDS preference, in the following, we only include main effects in the full and null models for comparisons.

4.5 NAE model (main effect of IDS preference)

Sys.setenv(R_FUTURE_RNG_ONMISUSE = "ignore") #this line of code is needed to get rid of warning message due to the "future" package. https://github.com/satijalab/seurat/issues/3622 

lmer.bf_full_1.nae <- brm(standardized.score.CDI.num ~ CDI.z_age_months + gender +
                        z_age_months + method + IDS_pref +
                        (1 | labid) + (1 | subid_unique),
                      data = data.nae, 
                      family = gaussian,
                      prior = prior_nae_1,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed=456)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
## 
## SAMPLING FOR MODEL '6361f04d0b2eda23c05b2e4be9a58962' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 7.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.75 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 2.19751 seconds (Warm-up)
## Chain 1:                2.05293 seconds (Sampling)
## Chain 1:                4.25045 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '6361f04d0b2eda23c05b2e4be9a58962' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 5.9e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.59 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 2.35878 seconds (Warm-up)
## Chain 2:                2.07516 seconds (Sampling)
## Chain 2:                4.43394 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '6361f04d0b2eda23c05b2e4be9a58962' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 6.4e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.64 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 2.25339 seconds (Warm-up)
## Chain 3:                1.3515 seconds (Sampling)
## Chain 3:                3.60489 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '6361f04d0b2eda23c05b2e4be9a58962' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 6.2e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.62 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 2.31966 seconds (Warm-up)
## Chain 4:                2.0556 seconds (Sampling)
## Chain 4:                4.37526 seconds (Total)
## Chain 4:
summary(lmer.bf_full_1.nae)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + (1 | labid) + (1 | subid_unique) 
##    Data: data.nae (Number of observations: 307) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 8) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.05      0.75     0.05     2.78 1.00     5477     4927
## 
## ~subid_unique (Number of levels: 211) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.87      0.66     0.03     2.44 1.00     6047     4451
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           40.24      2.10    36.10    44.27 1.00    11676     5899
## CDI.z_age_months     0.14      0.55    -0.92     1.20 1.00    16386     5805
## gender1             -1.98      1.61    -5.11     1.19 1.00    15997     5403
## z_age_months        -0.09      0.43    -0.93     0.75 1.00    14152     5721
## method1              0.19      1.28    -2.31     2.96 1.00    10776     4501
## method2             -0.05      1.79    -3.86     3.55 1.00     8843     3907
## IDS_pref            -0.62      2.16    -6.38     2.91 1.00     7741     3066
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    28.08      1.17    25.89    30.48 1.00    13999     5427
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_1.nae <- brm(standardized.score.CDI.num ~ CDI.z_age_months + gender +
                        z_age_months + method + 
                        (1 | labid) + (1 | subid_unique),
                      data = data.nae, 
                      family = gaussian,
                      prior = prior_nae_null_1,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed =123)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
## 
## SAMPLING FOR MODEL '4bb864c3fbd97441232e09115f930fe8' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 8.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.81 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 2.40827 seconds (Warm-up)
## Chain 1:                1.76915 seconds (Sampling)
## Chain 1:                4.17742 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '4bb864c3fbd97441232e09115f930fe8' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.00012 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.2 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 2.0687 seconds (Warm-up)
## Chain 2:                2.04309 seconds (Sampling)
## Chain 2:                4.11179 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '4bb864c3fbd97441232e09115f930fe8' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000298 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 2.98 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 2.14399 seconds (Warm-up)
## Chain 3:                2.02593 seconds (Sampling)
## Chain 3:                4.16992 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '4bb864c3fbd97441232e09115f930fe8' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 7.9e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.79 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 2.128 seconds (Warm-up)
## Chain 4:                2.0381 seconds (Sampling)
## Chain 4:                4.1661 seconds (Total)
## Chain 4:
summary(lmer.bf_null_1.nae)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months + method + (1 | labid) + (1 | subid_unique) 
##    Data: data.nae (Number of observations: 307) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 8) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.05      0.74     0.05     2.74 1.00     6284     4654
## 
## ~subid_unique (Number of levels: 211) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.86      0.66     0.03     2.42 1.00     5865     4466
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           40.14      2.11    35.93    44.29 1.00    13214     5661
## CDI.z_age_months     0.14      0.55    -0.93     1.22 1.00    16648     5387
## gender1             -1.97      1.65    -5.17     1.24 1.00    17327     5420
## z_age_months        -0.11      0.42    -0.95     0.71 1.00    16529     5611
## method1              0.22      1.22    -2.15     2.89 1.00    10974     4928
## method2             -0.09      1.81    -4.07     3.51 1.00     9396     3756
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    28.10      1.15    25.97    30.42 1.00    17149     5479
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Calculate the marginal log likelihoods for each hypothese using Bridge sampling

marloglik_full_1.nae <- bridge_sampler(lmer.bf_full_1.nae, silent = T)
marloglik_null_1.nae <- bridge_sampler(lmer.bf_null_1.nae, silent = T)

Bayes factor

BF_full_1.nae <- bayes_factor(marloglik_full_1.nae, marloglik_null_1.nae)
BF_full_1.nae
## Estimated Bayes factor in favor of x1 over x2: 0.63608

4.6 NAE model (interaction between IDS test age and IDS preference)

lmer.bf_full_2.nae <- brm(standardized.score.CDI.num ~ CDI.z_age_months + gender +
                        z_age_months + method + IDS_pref +
                        IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
                        (1 | labid) + (1 | subid_unique),
                      data = data.nae, 
                      family = gaussian,
                      prior = prior_nae_2,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed=890)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
## 
## SAMPLING FOR MODEL 'f3788e4416fd7ee1c5e7b74b10adc8b6' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000157 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.57 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 2.48097 seconds (Warm-up)
## Chain 1:                2.15046 seconds (Sampling)
## Chain 1:                4.63143 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'f3788e4416fd7ee1c5e7b74b10adc8b6' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000325 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 3.25 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 2.38087 seconds (Warm-up)
## Chain 2:                2.13843 seconds (Sampling)
## Chain 2:                4.5193 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'f3788e4416fd7ee1c5e7b74b10adc8b6' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000124 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.24 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 2.27497 seconds (Warm-up)
## Chain 3:                2.16732 seconds (Sampling)
## Chain 3:                4.44229 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'f3788e4416fd7ee1c5e7b74b10adc8b6' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 8.8e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.88 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 2.49946 seconds (Warm-up)
## Chain 4:                2.14814 seconds (Sampling)
## Chain 4:                4.64759 seconds (Total)
## Chain 4:
summary(lmer.bf_full_2.nae)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 | labid) + (1 | subid_unique) 
##    Data: data.nae (Number of observations: 307) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 8) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.03      0.75     0.04     2.78 1.00     5677     3586
## 
## ~subid_unique (Number of levels: 211) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.89      0.66     0.04     2.46 1.00     7136     4836
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                    40.08      2.23    35.69    44.40 1.00     8169
## CDI.z_age_months              0.08      0.58    -1.09     1.23 1.00    18071
## gender1                      -1.91      1.61    -5.03     1.27 1.00    19483
## z_age_months                 -0.18      0.44    -1.07     0.68 1.00    18190
## method1                       0.23      1.38    -2.36     3.19 1.00     7189
## method2                      -0.16      2.16    -4.71     3.72 1.00     5468
## IDS_pref                     -0.91      2.51    -7.88     2.79 1.00     7484
## method1:IDS_pref             -0.01      2.04    -4.54     4.43 1.00    11545
## method2:IDS_pref              0.35      2.54    -4.30     6.75 1.00     7882
## CDI.z_age_months:IDS_pref     0.45      1.07    -1.44     2.90 1.00    13628
## z_age_months:IDS_pref         0.81      1.15    -1.11     3.49 1.00    13396
##                           Tail_ESS
## Intercept                     4067
## CDI.z_age_months              5831
## gender1                       5281
## z_age_months                  5450
## method1                       2950
## method2                       2415
## IDS_pref                      3357
## method1:IDS_pref              3582
## method2:IDS_pref              2938
## CDI.z_age_months:IDS_pref     5707
## z_age_months:IDS_pref         5809
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    28.09      1.17    25.92    30.45 1.00    19893     5820
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_2.nae <- brm(standardized.score.CDI.num ~ CDI.z_age_months + gender +
                        z_age_months + method + IDS_pref +
                        IDS_pref:method + IDS_pref:CDI.z_age_months + 
                        (1 | labid) + (1 | subid_unique),
                      data = data.nae, 
                      family = gaussian,
                      prior = prior_nae_null_2,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed =102)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
## 
## SAMPLING FOR MODEL '0e2a1080bab56c0e90a24e05388bf5a7' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000255 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.55 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 2.42628 seconds (Warm-up)
## Chain 1:                2.09945 seconds (Sampling)
## Chain 1:                4.52572 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '0e2a1080bab56c0e90a24e05388bf5a7' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000108 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.08 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 2.28607 seconds (Warm-up)
## Chain 2:                2.11654 seconds (Sampling)
## Chain 2:                4.40261 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '0e2a1080bab56c0e90a24e05388bf5a7' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000263 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 2.63 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 2.48995 seconds (Warm-up)
## Chain 3:                2.13717 seconds (Sampling)
## Chain 3:                4.62712 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '0e2a1080bab56c0e90a24e05388bf5a7' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000267 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 2.67 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 2.32855 seconds (Warm-up)
## Chain 4:                2.10664 seconds (Sampling)
## Chain 4:                4.43519 seconds (Total)
## Chain 4:
summary(lmer.bf_null_2.nae)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: standardized.score.CDI.num ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + (1 | labid) + (1 | subid_unique) 
##    Data: data.nae (Number of observations: 307) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 8) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.04      0.75     0.04     2.71 1.00     5908     4221
## 
## ~subid_unique (Number of levels: 211) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.88      0.65     0.04     2.41 1.00     6257     4689
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                    40.17      2.10    35.97    44.34 1.00    12196
## CDI.z_age_months              0.08      0.57    -1.04     1.21 1.00    17288
## gender1                      -1.97      1.60    -5.09     1.16 1.00    20629
## z_age_months                 -0.10      0.42    -0.95     0.75 1.00    17914
## method1                       0.21      1.29    -2.22     3.01 1.00    10045
## method2                      -0.10      1.90    -4.22     3.83 1.00     8727
## IDS_pref                     -0.61      2.16    -6.38     2.75 1.00     8840
## method1:IDS_pref             -0.07      2.00    -4.55     4.17 1.00    12060
## method2:IDS_pref              0.14      2.36    -4.77     5.66 1.00     8495
## CDI.z_age_months:IDS_pref     0.47      1.06    -1.44     2.86 1.00    15876
##                           Tail_ESS
## Intercept                     5688
## CDI.z_age_months              6153
## gender1                       6037
## z_age_months                  5889
## method1                       4185
## method2                       3640
## IDS_pref                      3328
## method1:IDS_pref              3278
## method2:IDS_pref              3061
## CDI.z_age_months:IDS_pref     5247
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    28.11      1.15    25.96    30.42 1.00    19020     6352
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Calculate the marginal log likelihoods for each hypothese using Bridge sampling

marloglik_full_2.nae <- bridge_sampler(lmer.bf_full_2.nae, silent = T)
marloglik_null_2.nae <- bridge_sampler(lmer.bf_null_2.nae, silent = T)

Bayes factor

BF_full_2.nae <- bayes_factor(marloglik_full_2.nae, marloglik_null_2.nae)
BF_full_2.nae
## Estimated Bayes factor in favor of x1 over x2: 0.98500

4.7 UK model (main effect of IDS preference)

lmer.bf_full_1.uk <- brm(vocab_nwords ~ CDI.z_age_months + gender +
                       z_age_months + method + IDS_pref +
                      (1 | labid) + (1 | subid_unique),
                     data = data.uk, 
                      family = gaussian,
                      prior = prior_uk_1,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed=213)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
## 
## SAMPLING FOR MODEL '379ff944c7b93801a64060c542798150' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 7e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.7 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 1.17577 seconds (Warm-up)
## Chain 1:                0.925538 seconds (Sampling)
## Chain 1:                2.10131 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '379ff944c7b93801a64060c542798150' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 7.9e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.79 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.896499 seconds (Warm-up)
## Chain 2:                0.956413 seconds (Sampling)
## Chain 2:                1.85291 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '379ff944c7b93801a64060c542798150' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 7e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.7 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 1.18867 seconds (Warm-up)
## Chain 3:                0.946536 seconds (Sampling)
## Chain 3:                2.1352 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '379ff944c7b93801a64060c542798150' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 6.1e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.61 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 1.05547 seconds (Warm-up)
## Chain 4:                1.53038 seconds (Sampling)
## Chain 4:                2.58585 seconds (Total)
## Chain 4:
summary(lmer.bf_full_1.uk)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + (1 | labid) + (1 | subid_unique) 
##    Data: data.uk (Number of observations: 119) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.81      0.61     0.04     2.28 1.00     5852     3538
## 
## ~subid_unique (Number of levels: 88) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.79      0.61     0.03     2.27 1.00     5496     3000
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          153.94      8.21   137.27   169.56 1.00     9730     5591
## CDI.z_age_months    27.35      2.94    21.38    32.85 1.00    12045     5528
## gender1              9.90      9.69    -8.52    28.64 1.00    10527     5291
## z_age_months        -1.03      1.63    -5.02     1.55 1.00     8414     4845
## method1             -0.69      3.17    -9.13     4.12 1.00     5119     1963
## IDS_pref             0.59      5.30    -7.24    13.01 1.00     4256     1457
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    91.43      6.38    80.20   105.07 1.00     9396     5966
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_1.uk <- brm(vocab_nwords ~ CDI.z_age_months + gender +
                       z_age_months + method + 
                      (1 | labid) + (1 | subid_unique),
                      data = data.uk, 
                      family = gaussian,
                      prior = prior_uk_null_1,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed =312)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
## 
## SAMPLING FOR MODEL '844f09eb03d0f7d774bc9ac00a910248' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.00023 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.3 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.917603 seconds (Warm-up)
## Chain 1:                0.853581 seconds (Sampling)
## Chain 1:                1.77118 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '844f09eb03d0f7d774bc9ac00a910248' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 4.1e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.41 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.769098 seconds (Warm-up)
## Chain 2:                0.95405 seconds (Sampling)
## Chain 2:                1.72315 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '844f09eb03d0f7d774bc9ac00a910248' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 6.2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.62 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.75625 seconds (Warm-up)
## Chain 3:                0.94152 seconds (Sampling)
## Chain 3:                1.69777 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '844f09eb03d0f7d774bc9ac00a910248' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000204 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 2.04 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.77893 seconds (Warm-up)
## Chain 4:                0.942255 seconds (Sampling)
## Chain 4:                1.72118 seconds (Total)
## Chain 4:
summary(lmer.bf_null_1.uk)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method + (1 | labid) + (1 | subid_unique) 
##    Data: data.uk (Number of observations: 119) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.81      0.60     0.03     2.22 1.00     6487     4278
## 
## ~subid_unique (Number of levels: 88) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.79      0.61     0.03     2.29 1.00     6551     3986
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          154.14      8.38   137.27   170.08 1.00    14105     5734
## CDI.z_age_months    27.43      2.93    21.44    33.06 1.00    16606     6113
## gender1              9.84      9.65    -8.35    28.43 1.00    15050     6205
## z_age_months        -1.02      1.64    -5.03     1.51 1.00     9120     4618
## method1             -0.71      3.19    -9.14     4.08 1.00     6670     2125
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    91.35      6.41    79.61   104.72 1.00    15410     6275
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Calculate the marginal log likelihoods for each hypothese using Bridge sampling

marloglik_full_1.uk <- bridge_sampler(lmer.bf_full_1.uk, silent = T)
marloglik_null_1.uk <- bridge_sampler(lmer.bf_null_1.uk, silent = T)

Bayes factor

BF_full_1.uk <- bayes_factor(marloglik_full_1.uk, marloglik_null_1.uk)
BF_full_1.uk
## Estimated Bayes factor in favor of x1 over x2: 0.94293

4.8 UK model (interaction between IDS test age and IDS preference)

lmer.bf_full_2.uk <- brm(vocab_nwords ~ CDI.z_age_months + gender +
                       z_age_months + method + IDS_pref +
                       IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
                      (1 | labid) + (1 | subid_unique),
                     data = data.uk, 
                      family = gaussian,
                      prior = prior_uk_2,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed=546)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
## 
## SAMPLING FOR MODEL '8936aa99d96ef130bfa122dd46b98c39' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 4.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.41 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 1.8092 seconds (Warm-up)
## Chain 1:                1.69476 seconds (Sampling)
## Chain 1:                3.50395 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '8936aa99d96ef130bfa122dd46b98c39' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.00025 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 2.5 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 1.69324 seconds (Warm-up)
## Chain 2:                0.992679 seconds (Sampling)
## Chain 2:                2.68592 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '8936aa99d96ef130bfa122dd46b98c39' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 6.7e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.67 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 1.689 seconds (Warm-up)
## Chain 3:                1.00774 seconds (Sampling)
## Chain 3:                2.69674 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '8936aa99d96ef130bfa122dd46b98c39' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 5.8e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.58 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 1.3114 seconds (Warm-up)
## Chain 4:                0.984855 seconds (Sampling)
## Chain 4:                2.29625 seconds (Total)
## Chain 4:
summary(lmer.bf_full_2.uk)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 | labid) + (1 | subid_unique) 
##    Data: data.uk (Number of observations: 119) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.81      0.61     0.04     2.26 1.00     5636     3152
## 
## ~subid_unique (Number of levels: 88) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.81      0.61     0.03     2.27 1.00     6346     3454
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                   154.11      8.17   137.45   169.50 1.00     8021
## CDI.z_age_months             27.41      2.93    21.50    33.06 1.00     8970
## gender1                      10.00      9.79    -8.47    29.05 1.00     8842
## z_age_months                 -1.04      1.66    -5.25     1.43 1.00     6597
## method1                      -0.69      3.10    -9.10     3.68 1.00     4440
## IDS_pref                      0.55      4.80    -6.98    11.94 1.00     3006
## method1:IDS_pref              0.01      4.24    -8.69     9.16 1.00     5044
## CDI.z_age_months:IDS_pref     0.05      2.29    -5.01     5.21 1.00     5910
## z_age_months:IDS_pref        -0.09      2.36    -5.13     4.70 1.00     5164
##                           Tail_ESS
## Intercept                     5564
## CDI.z_age_months              5616
## gender1                       5767
## z_age_months                  4177
## method1                       1560
## IDS_pref                      1117
## method1:IDS_pref              1968
## CDI.z_age_months:IDS_pref     2008
## z_age_months:IDS_pref         2415
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    91.50      6.54    79.89   105.11 1.00     9123     6176
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_2.uk <- brm(vocab_nwords ~ CDI.z_age_months + gender +
                       z_age_months + method + IDS_pref +
                       IDS_pref:method + IDS_pref:CDI.z_age_months + 
                      (1 | labid) + (1 | subid_unique),
                      data = data.uk, 
                      family = gaussian,
                      prior = prior_uk_null_2,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .95), 
                      seed =689)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
## 
## SAMPLING FOR MODEL '155a5bbbeeadf50ca6dd061d19360a9f' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 5.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.55 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 1.17176 seconds (Warm-up)
## Chain 1:                0.981056 seconds (Sampling)
## Chain 1:                2.15281 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '155a5bbbeeadf50ca6dd061d19360a9f' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 5.4e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.54 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 1.21728 seconds (Warm-up)
## Chain 2:                0.949116 seconds (Sampling)
## Chain 2:                2.1664 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '155a5bbbeeadf50ca6dd061d19360a9f' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 6.1e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.61 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 1.64771 seconds (Warm-up)
## Chain 3:                0.993318 seconds (Sampling)
## Chain 3:                2.64103 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '155a5bbbeeadf50ca6dd061d19360a9f' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000229 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 2.29 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 1.18755 seconds (Warm-up)
## Chain 4:                0.964777 seconds (Sampling)
## Chain 4:                2.15233 seconds (Total)
## Chain 4:
summary(lmer.bf_null_2.uk)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vocab_nwords ~ CDI.z_age_months + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + (1 | labid) + (1 | subid_unique) 
##    Data: data.uk (Number of observations: 119) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.80      0.60     0.03     2.26 1.00     7122     3819
## 
## ~subid_unique (Number of levels: 88) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.80      0.61     0.03     2.26 1.00     7561     3802
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                   154.24      8.33   137.43   170.45 1.00     9042
## CDI.z_age_months             27.39      2.98    21.38    33.10 1.00    10553
## gender1                      10.11      9.62    -8.08    28.88 1.00    11166
## z_age_months                 -1.04      1.68    -5.29     1.51 1.00     6649
## method1                      -0.78      3.41   -10.17     4.06 1.01     3586
## IDS_pref                      0.30      4.40    -6.96    10.83 1.00     4757
## method1:IDS_pref             -0.01      4.06    -8.00     8.28 1.00     4590
## CDI.z_age_months:IDS_pref     0.04      2.63    -5.12     5.61 1.01     3837
##                           Tail_ESS
## Intercept                     5701
## CDI.z_age_months              6000
## gender1                       6223
## z_age_months                  3744
## method1                       1188
## IDS_pref                      1472
## method1:IDS_pref              1860
## CDI.z_age_months:IDS_pref     1796
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    91.34      6.45    79.76   105.22 1.00    10269     6191
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Calculate the marginal log likelihoods for each hypothese using Bridge sampling

marloglik_full_2.uk <- bridge_sampler(lmer.bf_full_2.uk, silent = T)
marloglik_null_2.uk <- bridge_sampler(lmer.bf_null_2.uk, silent = T)

Bayes factor

BF_full_2.uk <- bayes_factor(marloglik_full_2.uk, marloglik_null_2.uk)
BF_full_2.uk
## Estimated Bayes factor in favor of x1 over x2: 0.92703

4.9 Combined model (interaction between lanaguage zone and IDS preference)

lmer.bf_full_3.combined <- brm(CDI.prop ~ CDI.z_age_months + language_zone + gender +
                           z_age_months + method + IDS_pref + IDS_pref:language_zone +
                           IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
                           (1 + CDI.z_age_months | labid) + (1 | subid_unique),
                         data = data.uk_nae,
                        family = gaussian,
                      prior = prior_combined_full_info,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .9999), 
                      seed=100)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
## 
## SAMPLING FOR MODEL 'dc68784226044421f8a0db87fb6a9761' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000201 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.01 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 18.3068 seconds (Warm-up)
## Chain 1:                7.8792 seconds (Sampling)
## Chain 1:                26.186 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'dc68784226044421f8a0db87fb6a9761' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000304 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 3.04 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 17.2098 seconds (Warm-up)
## Chain 2:                15.4955 seconds (Sampling)
## Chain 2:                32.7052 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'dc68784226044421f8a0db87fb6a9761' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000152 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.52 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 16.5475 seconds (Warm-up)
## Chain 3:                7.88142 seconds (Sampling)
## Chain 3:                24.429 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'dc68784226044421f8a0db87fb6a9761' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000151 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.51 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 16.5896 seconds (Warm-up)
## Chain 4:                15.5024 seconds (Sampling)
## Chain 4:                32.092 seconds (Total)
## Chain 4:
summary(lmer.bf_full_3.combined)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: CDI.prop ~ CDI.z_age_months + language_zone + gender + z_age_months + method + IDS_pref + IDS_pref:language_zone + IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 + CDI.z_age_months | labid) + (1 | subid_unique) 
##    Data: data.uk_nae (Number of observations: 424) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 13) 
##                                 Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                       0.01      0.01     0.00     0.02 1.00
## sd(CDI.z_age_months)                0.02      0.00     0.01     0.03 1.00
## cor(Intercept,CDI.z_age_months)     0.03      0.42    -0.76     0.79 1.00
##                                 Bulk_ESS Tail_ESS
## sd(Intercept)                       4550     3934
## sd(CDI.z_age_months)                6399     6415
## cor(Intercept,CDI.z_age_months)      973     2167
## 
## ~subid_unique (Number of levels: 302) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.01      0.01     0.00     0.03 1.00     4674     3846
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                     0.32      0.01     0.29     0.34 1.00    10446
## CDI.z_age_months              0.04      0.01     0.03     0.05 1.00     5095
## language_zone1                0.08      0.02     0.05     0.11 1.00     6294
## gender1                       0.03      0.01     0.01     0.05 1.00    12001
## z_age_months                 -0.00      0.00    -0.01     0.00 1.00    12834
## method1                      -0.00      0.01    -0.02     0.02 1.00     8983
## method2                       0.01      0.01    -0.02     0.04 1.00     6369
## IDS_pref                      0.00      0.02    -0.03     0.04 1.00    10120
## language_zone1:IDS_pref       0.01      0.02    -0.02     0.06 1.00     5807
## method1:IDS_pref             -0.00      0.02    -0.04     0.03 1.00     8979
## method2:IDS_pref             -0.00      0.02    -0.04     0.04 1.00     7079
## CDI.z_age_months:IDS_pref     0.00      0.01    -0.01     0.02 1.00    14388
## z_age_months:IDS_pref        -0.00      0.01    -0.02     0.01 1.00    12881
##                           Tail_ESS
## Intercept                     6105
## CDI.z_age_months              5274
## language_zone1                5396
## gender1                       5977
## z_age_months                  6187
## method1                       4958
## method2                       3343
## IDS_pref                      3866
## language_zone1:IDS_pref       2938
## method1:IDS_pref              3134
## method2:IDS_pref              3278
## CDI.z_age_months:IDS_pref     5442
## z_age_months:IDS_pref         5224
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.21      0.01     0.20     0.23 1.00    12956     5923
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
lmer.bf_null_3.combined <- brm(CDI.prop ~ CDI.z_age_months + language_zone + gender +
                           z_age_months + method + IDS_pref + 
                           IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months +
                           (1 + CDI.z_age_months | labid) + (1 | subid_unique),
                         data = data.uk_nae,
                      family = gaussian,
                      prior = prior_combined_null_info,
                      iter = 4000, 
                      warmup = 2000,
                      chains = 4,
                      save_pars = save_pars(all = T),
                      control = list(adapt_delta = .9999), 
                      seed =200)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
## 
## SAMPLING FOR MODEL 'f762272a91b9146804973af2bbae2eb4' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000515 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.15 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 18.0606 seconds (Warm-up)
## Chain 1:                15.7231 seconds (Sampling)
## Chain 1:                33.7837 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'f762272a91b9146804973af2bbae2eb4' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000183 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.83 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 17.0697 seconds (Warm-up)
## Chain 2:                15.3224 seconds (Sampling)
## Chain 2:                32.3921 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'f762272a91b9146804973af2bbae2eb4' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000333 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 3.33 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 18.7739 seconds (Warm-up)
## Chain 3:                8.02287 seconds (Sampling)
## Chain 3:                26.7968 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'f762272a91b9146804973af2bbae2eb4' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000525 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 5.25 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 18.542 seconds (Warm-up)
## Chain 4:                7.94283 seconds (Sampling)
## Chain 4:                26.4848 seconds (Total)
## Chain 4:
summary(lmer.bf_null_3.combined)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: CDI.prop ~ CDI.z_age_months + language_zone + gender + z_age_months + method + IDS_pref + IDS_pref:method + IDS_pref:CDI.z_age_months + IDS_pref:z_age_months + (1 + CDI.z_age_months | labid) + (1 | subid_unique) 
##    Data: data.uk_nae (Number of observations: 424) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~labid (Number of levels: 13) 
##                                 Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                       0.01      0.01     0.00     0.02 1.00
## sd(CDI.z_age_months)                0.02      0.00     0.01     0.03 1.00
## cor(Intercept,CDI.z_age_months)     0.01      0.43    -0.79     0.79 1.00
##                                 Bulk_ESS Tail_ESS
## sd(Intercept)                       5150     4137
## sd(CDI.z_age_months)                7566     6438
## cor(Intercept,CDI.z_age_months)     1200     2826
## 
## ~subid_unique (Number of levels: 302) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.01      0.01     0.00     0.03 1.00     4593     3446
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                     0.31      0.01     0.29     0.34 1.00    10940
## CDI.z_age_months              0.04      0.01     0.03     0.05 1.00     6321
## language_zone1                0.08      0.02     0.05     0.11 1.00     7368
## gender1                       0.03      0.01     0.01     0.05 1.00    15525
## z_age_months                 -0.00      0.00    -0.01     0.00 1.00    14549
## method1                      -0.00      0.01    -0.02     0.02 1.00    11072
## method2                       0.01      0.02    -0.02     0.04 1.00     7455
## IDS_pref                      0.00      0.02    -0.03     0.03 1.00    11181
## method1:IDS_pref             -0.00      0.02    -0.04     0.03 1.00     9539
## method2:IDS_pref              0.00      0.02    -0.03     0.04 1.00     9235
## CDI.z_age_months:IDS_pref     0.00      0.01    -0.01     0.02 1.00    15202
## z_age_months:IDS_pref        -0.00      0.01    -0.02     0.01 1.00    12744
##                           Tail_ESS
## Intercept                     6379
## CDI.z_age_months              5664
## language_zone1                5211
## gender1                       5569
## z_age_months                  6425
## method1                       5420
## method2                       3905
## IDS_pref                      4151
## method1:IDS_pref              3629
## method2:IDS_pref              3787
## CDI.z_age_months:IDS_pref     4821
## z_age_months:IDS_pref         5179
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.21      0.01     0.20     0.23 1.00    13563     6219
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Calculate the marginal log likelihoods for each hypothese using Bridge sampling

marloglik_full_3.combined <- bridge_sampler(lmer.bf_full_3.combined, silent = T)
marloglik_null_3.combined <- bridge_sampler(lmer.bf_null_3.combined, silent = T)

Bayes factor

BF_full_3.combined <- bayes_factor(marloglik_full_3.combined, marloglik_null_3.combined)
BF_full_3.combined
## Estimated Bayes factor in favor of x1 over x2: 0.06768

```